Sistemas y Señales Biomédicos

Ingeniería Biomédica

Ph.D. Pablo Eduardo Caicedo Rodríguez

2025-05-05

Sistemas y Señales Biomedicos - SYSB

Digital Filters

Why the Z-Transform?

  • The Fourier Transform assumes signals are stable and well-behaved
  • But some biosignals or systems may not be absolutely summable
  • The Z-Transform generalizes the Fourier Transform
  • Useful for analyzing discrete-time systems, especially when stability and causality matter

Definition

Let \(x[n]\) be a discrete-time signal.

The Z-Transform is defined as:

\[X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}\]

Where: - \(z \in \mathbb{C}\) is a complex variable - \(z = re^{j\omega}\)

Region of Convergence (ROC)

  • The Z-Transform converges only for certain values of \(z\)
  • The set of \(z\) for which the series converges is the ROC
  • ROC is critical for system stability and causality

Causal Signals
ROC is outside outermost pole

Anti-Causal Signals
ROC is inside innermost pole

Z-Plane Representation

  • Poles: values of \(z\) where \(X(z) \to \infty\)
  • Zeros: values where \(X(z) = 0\)
  • Visualization of poles and zeros helps in understanding system behavior

\[H(z) = 1.00 \cdot \frac{(z - 0.50)}{(z - 0.90)}\]

(-1.5, 1.5)
(-1.5, 1.5)

Relationship with Fourier Transform

If the ROC includes the unit circle, \(|z| = 1\), then:

\[X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}\]

So the Fourier Transform is a special case of the Z-Transform.

Example Signal

Let:

\[x[n] = a^n u[n]\]

Where: - \(a \in \mathbb{R}\) - \(u[n]\) is the unit step function (0 for \(n<0\), 1 for \(n\geq 0\))

Step 1: Z-Transform

Apply the definition:

\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \sum_{n=0}^{\infty} (az^{-1})^n\]

This is a geometric series:

\[X(z) = \frac{1}{1 - az^{-1}} = \frac{z}{z - a}\]

What is a Geometric Series?

A geometric series is a sum of terms where each term is multiplied by the same constant:

\[S = \sum_{n=0}^{\infty} r^n = 1 + r + r^2 + r^3 + \cdots\]

The value of \(r\) determines whether this sum converges (has a finite limit) or diverges.

Goal

Understand why this series:

\[\sum_{n=0}^{\infty} r^n\]

converges if and only if:

\[\boxed{|r| < 1}\]

Partial Sums

Let’s consider the sum up to the \(N\)-th term:

\[S_N = \sum_{n=0}^{N} r^n = 1 + r + r^2 + \cdots + r^N\]

This has a closed-form expression:

\[S_N = \frac{1 - r^{N+1}}{1 - r} \quad \text{for } r \neq 1\]

Taking the Limit

To find the sum of the infinite series, take the limit as \(N \to \infty\):

\[S = \lim_{N \to \infty} S_N = \lim_{N \to \infty} \frac{1 - r^{N+1}}{1 - r}\]

Case 1: \(|r| < 1\)

If \(|r| < 1\), then:

\[r^{N+1} \to 0 \quad \text{as } N \to \infty\]

So the sum becomes:

\[S = \frac{1}{1 - r}\]

✅ The geometric series converges.

Case 2: \(|r| \geq 1\)

  • If \(r = 1\), then \(S_N = N + 1 \to \infty\)
  • If \(r = -1\), the sum oscillates: \(1 - 1 + 1 - 1 + \cdots\)
  • If \(|r| > 1\), then \(r^{N+1} \to \infty\)

❌ In all cases: the series diverges

Intuition

  • When \(|r| < 1\), each term \(r^n\) gets smaller and smaller
  • Their total sum settles to a finite number
  • When \(|r| \geq 1\), the terms don’t vanish — the sum keeps growing or oscillating

Example: \(r = 0.5\)

\[S = 1 + 0.5 + 0.25 + 0.125 + \cdots = \frac{1}{1 - 0.5} = 2\]

Every term adds less. The sum “flattens out.”

Why This Matters

The Z-Transform often gives us geometric series like:

\[\sum_{n=0}^{\infty} (az^{-1})^n\]

This converges only if:

\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]

So, understanding convergence of geometric series = understanding ROC in Z-transforms.

Summary

Condition Behavior Result
\(|r| < 1\) Terms shrink Series converges
\(|r| \geq 1\) Terms grow or oscillate Diverges

\[\sum_{n=0}^{\infty} r^n = \frac{1}{1 - r} \quad \text{if } |r| < 1\]

Step 2: Region of Convergence (ROC)

For convergence of the geometric series:

\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]

Therefore, the ROC is:

\[\boxed{|z| > |a|}\]

Interpretation

  • Causal signal (defined for \(n \geq 0\))
  • ROC is outside the outermost pole
  • Stable system only if ROC includes the unit circle \(|z| = 1\)

Example: \(a = 0.5\)

\[x[n] = (0.5)^n u[n]\]

Z-Transform:

\[X(z) = \sum_{n=0}^{\infty} (0.5)^n z^{-n} = \frac{z}{z - 0.5}\]

Region of Convergence:

\[\boxed{|z| > 0.5}\]

Summary

Signal Z-Transform ROC
\(x[n] = a^n u[n]\) \(\frac{z}{z - a}\) \(|z| > |a|\)
Example: \(a = 0.5\) \(\frac{z}{z - 0.5}\) \(|z| > 0.5\)

Properties of the Z-Transform

  • Linearity: \(a x[n] + b y[n] \to aX(z) + bY(z)\)
  • Time shifting: \(x[n - k] \to z^{-k} X(z)\)
  • Scaling in the z-domain: \(a^n x[n] \to X(z/a)\)
  • Convolution: \(x[n] * h[n] \to X(z)H(z)\)

Example

Let \(x[n] = a^n u[n]\), where \(|a| < 1\)

\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \frac{1}{1 - az^{-1}}, \quad \text{ROC: } |z| > |a|\]

Difference Equations in DSP

A difference equation relates input and output values at different time steps.

\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]

Common in: - Digital filters (FIR, IIR) - Signal models in ECG, EEG analysis - Implementation in real-time biosignal systems

Z-Transform of Time-Shifted Terms

The Z-Transform turns time shifts into powers of \(z^{-1}\):

Time Domain Z-Domain
\(x[n]\) \(X(z)\)
\(x[n-k]\) \(z^{-k} X(z)\)
\(y[n-k]\) \(z^{-k} Y(z)\)

Step 1: Apply Z-Transform

Given:

\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]

Apply \(\mathcal{Z} \{ \cdot \}\):

\[Y(z) - a_1 z^{-1} Y(z) - a_2 z^{-2} Y(z) = b_0 X(z) + b_1 z^{-1} X(z)\]

Step 2: Factor and Solve for \(H(z)\)

Group:

\[Y(z)(1 - a_1 z^{-1} - a_2 z^{-2}) = X(z)(b_0 + b_1 z^{-1})\]

Divide both sides:

\[H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1}}{1 - a_1 z^{-1} - a_2 z^{-2}}\]

Example

Given:

\[y[n] - 0.9 y[n-1] = x[n] - 0.5 x[n-1]\]

Z-Transform:

\[Y(z)(1 - 0.9 z^{-1}) = X(z)(1 - 0.5 z^{-1})\]

Transfer Function:

\[H(z) = \frac{1 - 0.5 z^{-1}}{1 - 0.9 z^{-1}}\]

Poles and Zeros

Let’s analyze \(H(z)\):

  • Zeros: Roots of the numerator \(\Rightarrow z = 0.5\)
  • Poles: Roots of the denominator \(\Rightarrow z = 0.9\)

Pole-Zero Plot
Visualizes system behavior
Check for: - Stability (poles inside unit circle) - Frequency shaping

Practice

Convert this equation:

\[y[n] = 0.6 y[n-1] + x[n] + x[n-1]\]

Find: - \(H(z)\) - Poles and zeros - Plot them in the Z-plane

Application in Biosignal Processing

  • Analysis of digital filters for ECG, EEG, etc.
  • Design of stable and causal filtering systems
  • Useful in difference equation modeling of biosignals

Summary

  • Z-Transform is a powerful tool for analyzing discrete systems
  • Provides insight into stability, causality, and system behavior
  • A generalization of the Fourier Transform
  • Crucial in digital signal processing of biosignals
  • Z-Transform converts difference equations into algebraic expressions
  • Transfer function \(H(z)\) tells us how the system responds to inputs
  • Key for digital filter design in biosignal processing

Next Steps

  • Practice Z-Transform computations
  • Pole-zero plotting exercises
  • Application to real biosignal filtering problems

Filter Design

Notch Filter

A common issue in biosignal processing is removing power‐line interference (50/60 Hz) from, for example, EEG or ECG signals. A simple digital filter to eliminate 60 Hz interference (assuming a sampling frequency \(f_s = 5000\) Hz) is to place complex‐conjugate zeros at

\[ z = e^{\pm j 2\pi\frac{60}{5000}}. \]

The resulting transfer function can be written as

\[ H(z) = 1 \;-\; 2\cos\!\Bigl(2\pi\frac{60}{5000}\Bigr)\,z^{-1} \;+\; z^{-2}. \]

This \(H(z)\) has zeros at \(e^{\pm j2\pi(60/5000)}\) that precisely cancel the 60 Hz component, thereby implementing a notch filter. Moreover, it is a second‐order FIR filter with symmetric coefficients, which grants it linear phase (important to avoid waveform distortion).

Filter Design

def design_filter(zeros=None, poles=None, gain=1.0):
    """
    Diseña un filtro digital a partir de ceros y/o polos y una ganancia.

    Parámetros:
    - zeros: lista de ceros (raíces del numerador), o None para no incluir
    - poles: lista de polos (raíces del denominador), o None para no incluir
    - gain: ganancia escalar del filtro

    Devuelve:
    - b: coeficientes del numerador
    - a: coeficientes del denominador
    """
    # Si no se pasan ceros, asumimos un FIR trivial (b = [gain])
    if zeros:
        b = gain * np.poly(zeros)
    else:
        b = np.array([gain], dtype=float)

    # Si no se pasan polos, asumimos sistema FIR (a = [1])
    if poles:
        a = np.poly(poles)
    else:
        a = np.array([1.0], dtype=float)

    return b, a

Filter Design

def gaussian(x, mu, sigma, A):
    """
    Genera una función gaussiana.

    Parámetros:
    - x: array de tiempos
    - mu: posición central de la gaussiana
    - sigma: desviación estándar
    - A: amplitud
    """
    return A * np.exp(-((x - mu) ** 2) / (2 * sigma**2))

Filter Design

def simulate_ecg(duration=10.0, fs=500, heart_rate=60):
    """
    Simula un ECG sintético basado en la superposición de ondas gaussianas.

    Parámetros:
    - duration: duración de la señal en segundos
    - fs: frecuencia de muestreo en Hz
    - heart_rate: frecuencia cardiaca en latidos por minuto

    Devuelve:
    - t: vector de tiempos
    - ecg: señal simulada de ECG en mV
    """
    dt = 1 / fs
    t = np.arange(0, duration, dt)
    rr = 60 / heart_rate  # intervalo RR en segundos

    # Inicializar señal
    ecg = np.zeros_like(t)

    # Parámetros de las ondas (posiciones relativas en segundos)
    # P wave
    p_amp, p_dur, p_delay = 0.25, 0.09, 0.16
    # Q wave
    q_amp, q_dur, q_delay = -0.05, 0.066, 0.166
    # R wave
    r_amp, r_dur, r_delay = 1.6, 0.1, 0.166
    # S wave
    s_amp, s_dur, s_delay = -0.25, 0.066, 0.19
    # T wave
    t_amp, t_dur, t_delay = 0.35, 0.142, 0.36

    # Generar cada latido
    for beat_start in np.arange(0, duration, rr):
        mask = (t >= beat_start) & (t < beat_start + rr)
        tb = t[mask] - beat_start
        ecg[mask] += gaussian(tb, p_delay, p_dur / 2, p_amp)
        ecg[mask] += gaussian(tb, q_delay, q_dur / 2, q_amp)
        ecg[mask] += gaussian(tb, r_delay, r_dur / 2, r_amp)
        ecg[mask] += gaussian(tb, s_delay, s_dur / 2, s_amp)
        ecg[mask] += gaussian(tb, t_delay, t_dur / 2, t_amp)

    return t, ecg+0.1*np.cos(2*np.pi*60*t)

Filter Design

# Parámetros de simulación
DURATION = 10    # segundos
FS = 500         # Hz
HR = 70          # latidos por minuto

# Generar señal
t, ecg_signal = simulate_ecg(duration=DURATION, fs=FS, heart_rate=HR)

# Graficar resultado
plt.figure(figsize=(12, 4))
plt.plot(t, ecg_signal, linewidth=1)
plt.title(f'Señal de ECG sintética ({HR} bpm)')
plt.xlabel('Tiempo (s)')
plt.ylabel('Amplitud (mV)')
plt.grid(True)
plt.tight_layout()
plt.show()

Filter Design

n = len(ecg_signal)
yf = np.fft.fft(ecg_signal)
xf = np.fft.fftfreq(n, 1/fs)[:n//2]
plt.plot(xf, 2.0/n * np.abs(yf[0:n//2]))
plt.grid()
plt.show()

Filter Design

fc = 60
b,a =design_filter(zeros=[np.exp(1j*2*np.pi*fc/fs), np.exp(-1j*2*np.pi*fc/fs)])
print(a)
print(b)

w, h = sig.freqz(a, b, worN=8000)
plt.plot(w/np.pi*fs/2, 20*np.log10(abs(h)))
plt.grid()
plt.xlabel('Frequency (Hz)')

Filter Design

FIR filters

FIR filters (Finite Impulse Response) are widely used in biomedical processing because they can be designed to have linear phase response, avoiding phase distortion in the filtered signal (which is useful for preserving the morphology of ECG waves, for example).

Filter Design

Filter Design Process

  1. Defining the desired ideal frequency response \(H_d(e^{j\omega})\).

  2. Obtaining the ideal impulse response \(h_d[n]\) as the inverse Fourier transform of \(H_d\).

  3. Truncating \(h_d[n]\) (which is usually infinite or very long) with a window function \(w[n]\) to obtain a realizable FIR filter

    \[ h[n] = h_d[n]\,w[n]. \]

Filter Design

You obtain \(h_d[n]\) as the inverse discrete–time Fourier transform of your ideal frequency response \(H_d(e^{j\omega})\). For a low-pass filter with cutoff \(\omega_c\),

\[ H_d(e^{j\omega}) = \begin{cases} 1, & |\omega|\le\omega_c,\\ 0, & \text{otherwise}. \end{cases} \]

By definition of the inverse DTFT,

\[ h_d[n] \;=\; \frac{1}{2\pi}\int_{-\pi}^{\pi} H_d(e^{j\omega})\,e^{j\omega n}\,d\omega \;=\;\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} e^{j\omega n}\,d\omega. \]

Filter Design

Carry out the integral in two cases:

  1. For \(n\neq 0\):

    \[ h_d[n] =\frac{1}{2\pi}\,\frac{e^{j\omega n}}{j\,n}\Biggr|_{-\omega_c}^{\omega_c} =\frac{1}{2\pi}\,\frac{e^{j\omega_c n}-e^{-j\omega_c n}}{j\,n} =\frac{\sin(\omega_c\,n)}{\pi\,n}. \]

  2. For \(n=0\):

    \[ h_d[0] =\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} 1\,d\omega =\frac{2\,\omega_c}{2\pi} =\frac{\omega_c}{\pi}. \]

Putting both together,

\[ h_d[n] =\begin{cases} \dfrac{\sin(\omega_c\,n)}{\pi\,n}, & n\neq 0,\\[1em] \dfrac{\omega_c}{\pi}, & n=0. \end{cases} \]

Relating to sampling frequency

If your cutoff is specified in Hz, \(f_c\), and sampling rate is \(f_s\), then

\[ \omega_c = 2\pi\,\frac{f_c}{f_s}, \]

so you can write

\[ h_d[n] =\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{\pi\,n} =\;2\;\frac{f_c}{f_s}\;\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{2\pi \frac{f_c}{f_s}\,n} =2\frac{f_c}{f_s}\,\mathrm{sinc}\!\Bigl(2\frac{f_c}{f_s}\,n\Bigr), \]

where we define the normalized sinc as \(\mathrm{sinc}(x)=\frac{\sin(\pi x)}{\pi x}\).

Key takeaway

  • \(h_d[n]\) is exactly the inverse‐DTFT of the ideal (“brick‐wall”) frequency specification.
  • It produces a sinc-shaped impulse response of infinite length.
  • Truncation (with a window) makes it realizable as a finite-length FIR.

Filter Design

The typical characteristics of classic windows are:

  • Rectangular: narrowest main lobe (width ≈ 4π/M rad) but highest sidelobes (first sidelobe ≈ –13 dB, stop-band attenuation ≈ 21 dB). It gives the steepest transition for a given filter order, at the expense of poorer stop-band rejection.

  • Bartlett (triangular): somewhat wider main lobe (≈ 8π/M), sidelobes ≈ –25 dB.

  • Hann: main lobe ≈ 8π/M, sidelobes ≈ –31 dB (better rejection than rectangular but smoother transitions).

  • Hamming: main lobe ≈ 8π/M, sidelobes ≈ –41 dB (minimum stop-band attenuation ≈ 53 dB). Very popular for its good compromise between transition width and stop-band attenuation.

  • Blackman: wider main lobe (≈ 12π/M) but very low sidelobes (≈ –57 dB, attenuation ≈ 74 dB).

  • Kaiser: allows selection of a parameter β to control sidelobe attenuation, offering flexibility. Approximately, to achieve A dB of attenuation,

    \[ \beta \approx 0.1102\,(A - 8.7)\quad(\text{for }A>50), \]

    and the normalized transition width Δω relates to the order M and β by

    \[ M \approx \frac{A - 8}{2.285\,\Delta\omega} \]

These formulas stem from Kaiser’s approximations and help in sizing the filter.

Windows Forms

(-100.0, 5.0)

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import (
    boxcar,       # Rectangular
    bartlett,    # Triangular
    hann,        # Hann
    hamming,     # Hamming
    blackman,    # Blackman
    kaiser       # Kaiser
)

# Parameters
M = 64               # window length
beta = 8.6           # Kaiser parameter for moderate sidelobe attenuation
nfft = 512           # FFT length for frequency response
fs = 1.0             # normalized sampling rate

# Generate windows
windows = {
    'Rectangular': boxcar(M),
    'Bartlett':    bartlett(M),
    'Hann':        hann(M),
    'Hamming':     hamming(M),
    'Blackman':    blackman(M),
    f'Kaiser (β={beta})': kaiser(M, beta)
}

# Time-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
    plt.plot(np.arange(M), w, label=name)
plt.title('Window Functions — Time Domain')
plt.xlabel('Sample index n')
plt.ylabel('Amplitude')
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()

# Frequency-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
    # Compute normalized frequency response
    W = np.fft.fft(w, nfft)
    W_mag = 20 * np.log10(np.abs(W) / np.max(np.abs(W)))
    freqs = np.fft.fftfreq(nfft, d=1/fs)
    # Only plot 0 ≤ f ≤ 0.5 (normalized Nyquist)
    idx = np.logical_and(freqs >= 0, freqs <= 0.5)
    plt.plot(freqs[idx], W_mag[idx], label=name)
plt.title('Window Functions — Frequency Response')
plt.xlabel('Normalized Frequency (cycles/sample)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-100, 5)
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()

plt.show()

Filter Design

N = 101
fc = 30

# 2. Compute the ideal impulse response h_d[n] of a low-pass filter
n = np.arange(N)
M = (N - 1) / 2
# Using the normalized sinc: sinc(x) = sin(pi*x)/(pi*x)
h_d = (2 * fc / FS) * np.sinc(2 * fc * (n - M) / FS)

# 3. Choose a window w[n] (here: Hamming) and truncate
w = np.hamming(N)
h = h_d * w

# 4. Normalize to ensure unity gain at DC
h /= np.sum(h)

# 6. Filter it (only allowed library call)
y = sig.lfilter(h, 1.0, ecg_signal)

# 7. Plot input vs. output
plt.figure(figsize=(8, 4))
plt.plot(t, ecg_signal, label='Original signal')
plt.plot(t, y, label='Filtered signal', linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Low-pass FIR (window method) – Cutoff 100 Hz')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

Introduction: IIR Filter Design Strategy

  • Goal: Design digital Infinite Impulse Response (IIR) filters.
  • Approach: Leverage well-established theory and techniques from analog filter design.
  • Method:
    1. Design a suitable analog filter prototype \(H(s)\).
    2. Transform this analog design into a digital filter \(H(z)\).

Step 0: The Analog Domain - Laplace Transform

  • Context: Continuous-time signals and Linear Time-Invariant (LTI) systems.
  • Tool: The Laplace Transform converts differential equations (time domain) to algebraic equations (s-domain).
  • Transfer Function \(H(s)\): Represents the system in the s-domain. \[H(s) = \frac{Y(s)}{X(s)}\] where \(s = \sigma + j\Omega\) is the complex frequency (\(\Omega\) = analog angular frequency).

Step 0: The Analog Domain - s-Plane Analysis

  • s-Plane: Complex plane where \(H(s)\) is analyzed.
  • Poles & Zeros: Roots of the denominator & numerator of \(H(s)\), respectively. Their locations determine filter behavior.
  • Frequency Response: Determined by \(H(j\Omega)\) (evaluating \(H(s)\) on the imaginary axis).
  • Stability: For a causal system to be stable, all poles must be in the Left-Half Plane (LHP), i.e., \(\text{Re}\{s\} < 0\).

Step 0: The Analog Domain - Filter Prototypes

  • Standardized, well-understood analog filter approximations:
    • Butterworth: Maximally flat passband.
    • Chebyshev Type I: Equiripple passband, monotonic stopband.
    • Chebyshev Type II: Monotonic passband, equiripple stopband.
    • Elliptic (Cauer): Equiripple in both passband and stopband (sharpest transition for a given order).
  • Design usually starts with a normalized low-pass prototype.

The Bridge: Analog-to-Digital Mapping

  • Core Task: Convert the analog filter \(H(s)\) (s-plane) into a digital filter \(H(z)\) (z-plane).
  • Need a mapping that relates the continuous frequency \(\Omega\) to the discrete frequency \(\omega\).
  • Crucially, the stable region (LHP in s-plane) should map to the stable region (inside the unit circle in z-plane).

Step 1: Digital Filter Specifications

  • Define the desired characteristics of the digital filter:
    • Critical Frequencies:
      • Passband edge: \(\omega_p\)
      • Stopband edge: \(\omega_s\)
      • (Normalized digital angular frequencies, e.g., \(0 \le \omega \le \pi\))
    • Tolerances:
      • Max. passband ripple/attenuation: \(A_p\) (dB) or \(\delta_p\)
      • Min. stopband attenuation: \(A_s\) (dB) or \(\delta_s\)

Step 2: Transformation Methods - Overview

  • Mathematical mappings from the s-plane to the z-plane.
  • Two primary methods:
    1. Impulse Invariance
    2. Bilinear Transformation

Step 2: Method 1 - Impulse Invariance

  • Concept: Match the impulse response: \(h_{digital}[n] \approx T \cdot h_{analog}(nT)\).
  • Mapping: Maps s-plane pole \(s_k\) to z-plane pole \(z_k = e^{s_k T}\).
    • Uses partial fraction expansion of \(H(s)\).
  • Pros: Preserves impulse response shape.
  • Cons: Suffers from frequency aliasing if \(H(j\Omega)\) is not bandlimited below Nyquist frequency (\(F_s/2\)). Best for low-pass/band-pass filters.

Step 2: Method 2 - Bilinear Transformation

  • Concept: An algebraic substitution mapping the \(j\Omega\) axis onto the unit circle.
  • Mapping Formula: \[s = \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}\] (T = sampling period, often set to 2 for simplicity during derivation).
  • Pros:
    • No aliasing: Maps entire \(j\Omega\) axis to the unit circle (\(-\pi \le \omega \le \pi\)).
    • Preserves stability: Maps LHP (stable s-region) to inside the unit circle (stable z-region).
  • Cons: Introduces non-linear frequency warping.

Understanding Frequency Warping (Bilinear)

  • The relationship between analog (\(\Omega\)) and digital (\(\omega\)) frequencies is non-linear: \[\Omega = \frac{2}{T} \tan\left(\frac{\omega}{2}\right)\] or equivalently: \[\omega = 2 \arctan\left(\frac{\Omega T}{2}\right)\]
  • This compresses the infinite analog frequency axis onto the finite digital frequency range \([-\pi, \pi]\).
  • The mapping is non-uniform, most distorted near \(\omega = \pm \pi\).

Step 3: Frequency Pre-warping (Bilinear Only)

  • Problem: Frequency warping distorts the locations of \(\omega_p\) and \(\omega_s\).
  • Solution: Pre-warp the digital specifications (\(\omega_p, \omega_s\)) into corresponding analog frequencies (\(\Omega_p, \Omega_s\)) before designing the analog filter.
  • Formula: Use the warping relationship: \[\Omega_p = \frac{2}{T} \tan\left(\frac{\omega_p}{2}\right)\] \[\Omega_s = \frac{2}{T} \tan\left(\frac{\omega_s}{2}\right)\]
  • Use these \(\Omega_p, \Omega_s\) values for the analog design step.

Step 4: Analog Prototype Design Workflow

  1. Select Prototype: Choose Butterworth, Chebyshev, Elliptic based on specs (ripple, transition width) using the (pre-warped) \(\Omega_p, \Omega_s, A_p, A_s\).
  2. Determine Order (N) & Cutoff (\(\Omega_c\)): Calculate required analog filter order and cutoff frequency using prototype-specific formulas.
  3. Find Normalized LP \(H_{LP}(s)\): Obtain the transfer function for the normalized low-pass prototype.
  4. Frequency Transformation (Analog): If needed, transform \(H_{LP}(s)\) to target type (HP, BP, BS) using analog transformations.
    • Result: The final analog filter transfer function \(H_a(s)\).

Step 5: Apply Digital Transformation

  • Convert the designed analog filter \(H_a(s)\) to the digital domain \(H(z)\):
    • If Impulse Invariance: Use partial fractions and the \(s_k \rightarrow e^{s_k T}\) mapping.
    • If Bilinear Transformation: Substitute \(s = \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}\) into \(H_a(s)\) and simplify algebraically.
  • Result: The digital filter transfer function: \[H(z) = \frac{\sum_{k=0}^{M} b_k z^{-k}}{1 + \sum_{k=1}^{N} a_k z^{-k}}\]

Step 6: Realization (Implementation)

  • The transfer function \(H(z)\) corresponds directly to a difference equation: \[y[n] = \sum_{k=0}^{M} b_k x[n-k] - \sum_{k=1}^{N} a_k y[n-k]\]
  • \(x[n]\): Input sequence
  • \(y[n]\): Output sequence
  • \(b_k, a_k\): Filter coefficients derived from \(H(z)\).
  • The feedback term (sum involving \(y[n-k]\)) makes the impulse response potentially infinite (IIR).

Summary: IIR Design Process

  1. Define Digital Specs (\(\omega_p, \omega_s, A_p, A_s\)).
  2. Select Transformation (e.g., Bilinear).
  3. Pre-warp Frequencies (if Bilinear: \(\omega \rightarrow \Omega\)).
  4. Design Analog Filter \(H_a(s)\) (using prototypes & transformations).
  5. Transform to Digital \(H(z)\) (apply Bilinear or Impulse Invariance).
  6. Realize as Difference Equation for Implementation.

Introduction: IIR Filter Design Strategy

  • Goal: Design digital Infinite Impulse Response (IIR) filters.
  • Approach: Leverage well-established theory and techniques from analog filter design.
  • Method:
    1. Define digital specifications.
    2. (Conceptually) Design a suitable analog filter prototype \(H(s)\).
    3. (Conceptually) Transform \(H(s)\) into a digital filter \(H(z)\).
  • Practice (Code): Use libraries like scipy.signal that often encapsulate steps 2 & 3.

Step 0: The Analog Domain - Laplace Transform

  • Context: Continuous-time signals and Linear Time-Invariant (LTI) systems.
  • Tool: The Laplace Transform converts differential equations (time domain) to algebraic equations (s-domain).
  • Transfer Function \(H(s)\): Represents the system in the s-domain. \[H(s) = \frac{Y(s)}{X(s)}\] where \(s = \sigma + j\Omega\) is the complex frequency (\(\Omega\) = analog angular frequency).

Step 0: The Analog Domain - s-Plane Analysis

  • s-Plane: Complex plane where \(H(s)\) is analyzed.
  • Poles & Zeros: Roots of the denominator & numerator of \(H(s)\). Determine filter behavior.
  • Frequency Response: \(H(j\Omega)\).
  • Stability: All poles must be in the Left-Half Plane (LHP) (\(\text{Re}\{s\} < 0\)).

Step 0: The Analog Domain - Filter Prototypes

  • Standardized analog filter approximations:
    • Butterworth: Maximally flat passband.
    • Chebyshev I/II: Equiripple passband/stopband.
    • Elliptic (Cauer): Equiripple passband & stopband (sharpest transition).
  • These form the basis for many IIR digital filter designs.

The Bridge: Analog-to-Digital Mapping

  • Core Task: Convert \(H(s)\) (s-plane) \(\rightarrow\) \(H(z)\) (z-plane).
  • Relates continuous frequency \(\Omega\) to discrete frequency \(\omega\).
  • Stable LHP (s-plane) \(\rightarrow\) Stable interior of unit circle (z-plane).
  • Common Method: Bilinear Transformation (preserves stability, avoids aliasing, but warps frequency).

Step 1: Digital Filter Specifications

  • Define the desired characteristics of the digital filter:
    • Sampling Frequency: \(F_s\) (Hz)
    • Critical Frequencies (Digital):
      • Passband edge: \(f_p\) (Hz) \(\implies \omega_p = 2\pi f_p / F_s\) (radians/sample)
      • Stopband edge: \(f_s\) (Hz) \(\implies \omega_s = 2\pi f_s / F_s\) (radians/sample)
      • Note: Scipy often uses frequencies normalized to Nyquist: \(W_p = f_p / (F_s/2)\), \(W_s = f_s / (F_s/2)\)
    • Tolerances:
      • Max. passband ripple/attenuation: \(A_p\) or gpass (dB)
      • Min. stopband attenuation: \(A_s\) or gstop (dB)

Step 1: Example Specifications (Python)

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# Filter Specifications
fs = 10000  # Sampling frequency (Hz)
fp = 1000   # Passband edge frequency (Hz)
fs_stop = 1500 # Stopband edge frequency (Hz)
gpass = 1    # Passband ripple (dB) - max loss
gstop = 40   # Stopband attenuation (dB) - min attenuation

# Normalize frequencies to Nyquist frequency (fs/2)
wp = fp / (fs / 2)
ws = fs_stop / (fs / 2)

print(f"Sampling Frequency: {fs} Hz")
print(f"Normalized Passband Edge: {wp:.3f} (pi rad/sample)")
print(f"Normalized Stopband Edge: {ws:.3f} (pi rad/sample)")
print(f"Max Passband Loss: {gpass} dB")
print(f"Min Stopband Attenuation: {gstop} dB")

Steps 2-4: Finding Order & Designing (Conceptual)

  • Transformation: Bilinear Transform is implicitly used by many scipy.signal functions for IIR design.
  • Frequency Pre-warping: Handled internally by Scipy when given digital frequency specs.
  • Analog Design: Scipy calculates the required analog prototype order and parameters based on the (internally pre-warped) specs.
  • Goal: Find filter order \(N\) and digital cutoff frequency \(W_n\).

Steps 2-4: Finding Order & Wn (Python - Butterworth Example)

  • Use scipy.signal.buttord to find the minimum order and cutoff frequency for a digital Butterworth filter meeting the specs.
# --- Continued from previous code block ---

# Find the minimum order and cutoff frequency
N, Wn = signal.buttord(wp, ws, gpass, gstop, analog=False)

print(f"\nMinimum Filter Order (N): {N}")
# Wn is the digital cutoff freq (normalized) for butter function
print(f"Butterworth Natural Frequency (Wn): {Wn:.3f} (pi rad/sample)")
  • analog=False specifies we want parameters for a digital filter.
  • Wn is the digital cutoff frequency (scalar for LP/HP, tuple for BP/BS) required by the butter function. It’s often close to wp.

Step 5: Designing the Digital Filter (Python - Butterworth)

  • Use scipy.signal.butter (or cheby1, cheby2, ellip) to get the filter coefficients \(b\) (numerator) and \(a\) (denominator) of \(H(z)\).
# --- Continued from previous code block ---

# Design the digital Butterworth filter
# btype='low' for lowpass
b, a = signal.butter(N, Wn, btype='low', analog=False, output='ba')

print("\nFilter Coefficients:")
print(f"Numerator (b): {np.round(b, 5)}")
print(f"Denominator (a): {np.round(a, 5)}")
# H(z) = (b[0] + b[1]z^-1 + ...)/(a[0] + a[1]z^-1 + ...)
# Note: a[0] is usually 1 after normalization
  • output='ba' gives numerator/denominator coefficients. Other options: 'sos' (second-order sections - numerically better), 'zpk' (zeros, poles, gain).

Visualizing Frequency Response (Python)

  • Calculate the frequency response using scipy.signal.freqz.
  • Plot magnitude (in dB) and phase.
# --- Continued from previous code block ---

# Calculate frequency response
w, h = signal.freqz(b, a, worN=8000, fs=fs) # Use fs for physical freq axis

# Plot Magnitude Response
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(w, 20 * np.log10(abs(h)))
plt.title('Digital Butterworth Filter Frequency Response')
plt.ylabel('Amplitude [dB]')
plt.grid(True)
plt.axvline(fp, color='green', linestyle='--', label=f'fp={fp}Hz') # Passband edge
plt.axvline(fs_stop, color='red', linestyle='--', label=f'fs={fs_stop}Hz') # Stopband edge
plt.axhline(-gpass, color='green', linestyle=':', label=f'-{gpass}dB')
plt.axhline(-gstop, color='red', linestyle=':', label=f'-{gstop}dB')
plt.legend()
plt.xlim(0, fs/2)
plt.ylim(-60, 5) # Adjust ylim as needed

# Plot Phase Response
plt.subplot(2, 1, 2)
angles = np.unwrap(np.angle(h))
plt.plot(w, np.degrees(angles))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Angle [degrees]')
plt.grid(True)
plt.xlim(0, fs / 2)

plt.tight_layout()
plt.show()

Step 6: Realization & Filtering (Python Example)

  • The coefficients b, a define the difference equation.
  • Use scipy.signal.lfilter to apply the filter to a signal.
# --- Create a sample signal ---
T = 1.0 # seconds
n_samples = int(T * fs)
t = np.linspace(0, T, n_samples, endpoint=False)

# Signal: Sine wave inside passband + sine wave inside stopband
sig_pass = 0.8 * np.sin(2 * np.pi * fp/2 * t) # 500 Hz
sig_stop = 0.4 * np.sin(2 * np.pi * (fs_stop + 500) * t) # 2000 Hz
signal_in = sig_pass + sig_stop

# --- Apply the filter ---
filtered_signal = signal.lfilter(b, a, signal_in)

# --- Plot signals ---
plt.figure(figsize=(10, 5))
plt.plot(t, signal_in, label='Input Signal', alpha=0.7)
plt.plot(t, filtered_signal, label='Filtered Signal', linewidth=2)
plt.title('Filtering Example')
plt.xlabel('Time [sec]')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.ylim(-1.5, 1.5)
plt.show()

# Observe: The 500Hz component passes, the 2000Hz component is attenuated.

Summary: IIR Design Process (with Scipy)

  1. Define Digital Specs (\(f_p, f_s, gpass, gstop, F_s\)). Normalize frequencies (\(W_p, W_s\)).
  2. Select Filter Type (e.g., Butterworth) & find order/cutoff using buttord (digital).
  3. Design Filter to get coefficients (\(b, a\)) using butter (digital).
  4. Analyze Frequency Response using freqz.
  5. Implement/Apply Filter using lfilter (uses the difference equation implicitly).
  • Scipy handles the underlying analog prototype mapping, pre-warping, and bilinear transformation internally.